Method for calculating the electronic structure of correlated 
materials from a truly first-principles LDA+U scheme 

K. Karlsson^'^, F. Aryasetiawan^'^ and O. Jepsen^, 
^Department of Life Sciences, Hdgskolan i Skovde, 54128 Skovde, Sweden, 
"^Max Planck Institut fiir Festkorperforschung, D-705 06 Stuttgart, 
Germany, Ja.pa.n Science and Technology Agency, 
CREST, Kawaguchi 332-0012 Japan 
^Graduate School of Advanced Integration Science, 
Ghiba University, Ghiba, 263-8522 Japan, and 
Japan Science and Technology Agency, 
GREST, Kawaguchi 332-0012 Japan. 
(Dated:) 

Abstract 

We present a method for calculating the electronic structure of correlated materials based 
on a truly first-principles LDA+C/ scheme. Recently we suggested how to calculate U from 
first-principles, using a method which we named constrained RPA (cRPA). The input is simply 
the Kohn-Sham eigenfunctions and eigenvalues obtained within the LDA. In our proposed self- 
consistent LDA+f/ scheme, we calculate the LDA-I-C/ eigenfunctions and eigenvalues and use these 
to extract U . The updated U is then used in the next iteration to obtain a new set of eigenfunctions 
and eigenvalues and the iteration is continued until convergence is achieved. The most significant 
result is that our numerical approach is indeed stable: it is possible to find the effective exchange 
and correlation interaction matrix in a self- consistent way, resulting in a significant improvement 
over the LDA results, regarding both the bandgap in NiO and the /-band exchange spin-splitting 
in Gd, but some discrepancies still remain. 
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I. INTRODUCTION 



The interest for a fundamental understanding of strongly correlated systems has led to 
the development of a number of electronic structure methods. Among the most successful 
are the LDA+U approach proposed by Anisimov and coworkers [l| and the dynamical mean- 
field theory (DMFT) proposed by Georges and coworkers 2,1^. Very recently a new scheme, 
dubbed the LDA+Gutzwiller method .4|, for treating strong electron correlations was intro- 
duced. In all these methods the strong Coulomb onsite correlations for electrons residing 
in the localized orbitals are explicitly taken care of via a set of Hubbard like parameters 
or the Hubbard U. This is evidently unsatisfactory from the point of view of quantitative 
prediction of materials properties, since optical and magnetic excitations are of vital impor- 
tance in many technological applications such as solar cell design, optical memories, photo 
luminescent devices (semiconductor lasers and diodes) and photo chemical reactions. Often 
it has been shown that by adjusting the Hubbard U one can get results in good agreement 
with experiment but not for a good reason. Hereby lies the importance of determining U 
entirely from first-principles. 

Over the last two decades a number of methods for calculating the Hubbard U from 
first principles have been proposed. The pioneering work may be traced back to the paper 
by Gunnarsson et al. jsl who proposed to calculate U using the constrained LDA (cLDA) 
scheme. A few years ago, a new method for calculating the Hubbard [/, named the con- 
strained random-phase approximation (cRPA) method in analogy to the cLDA method, was 
proposed Q . The method allows for a systematic and precise determination of the Hubbard 
U entirely from first-principles from the knowledge of the bandstructure alone. The method 
was based on the intuitive idea that the Hubbard U can be viewed as a Coulomb interaction 
screened by the polarization of the whole system excluding the polarization arising from a 
set of bands which are treated in the Hubbard model. In other words, the Hubbard U when 
further screened by the electrons in the Hubbard model yields the screened interaction of the 
full system. This intuitive idea was recently shown to be rigorously correct and the cRPA is 
just an approximate way of calculating the screened interaction U within the random-phase 
approximation (RPA). 

By determining the Hubbard U from first-principles the cRPA method offers the pos- 
sibility of making methods based on the Hubbard U fully first-principles schemes. The 
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purpose of the present work is to develop a scheme for calculating the electronic structure 
of correlated materials based on a truly self- consistent first-principles LDA+f/ scheme. In 
conventional LDA+[/ scheme as it was originally proposed the Hubbard U is taken as an 
adjustable parameter which is fixed for a given calculation. In our proposed self-consistent 
LDA+f/ scheme, we calculate the LDA+t/ eigenfunctions and eigenvalues and use these to 
calculate U using the cRPA method. The new U is then used in the next iteration to obtain a 
new set of eigenfunctions and eigenvalues and the iteration is continued until convergence is 
achieved. Thus, U is no longer an arbitrarily adjusted parameter like in the original LDA+[/ 
scheme but rather it is determined self-consistently within the theoretical scheme. Our first 
target will be to calculate the electronic structure of the transition metal oxide series and as 
a test case we consider NiO, which is regarded as the epitome of Mott-Hubbard insulators. 
We are also aiming at obtaining a more satisfactory description of the electronic structure of 
the 4/ electron series which is highly problematic for the LDA. The path is then opened for 
more complex materials, such as magnetic semiconductors, for which no realistic methods 
are in existent at present. 

In the present paper we present some results for NiO and Gd, which we believe should 
provide us with a stringent test of the applicability of our method. The most important 
finding is that our numerical approach is indeed stable i.e it is possible to find U and J self- 
consistently. The bandgap in NiO and the spin-splitting of the /-bands in Gd are found to 
compare well with experiment using our self-consistent determined values of the correlation 
parameters. 

II. THEORY 



A. Constrained RPA 

We first give a short summary of the cRPA method presented in detail elsewhere js, 7| 
The fully screened Coulomb interaction is given by 



W=[l-vP]-\ (1) 
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where v is the bare Coulomb interaction and P is the non-interacting polarization given by 

occ unocc 
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where {ipijEi} are one-particle Bloch eigenfunctions and eigenvalues corresponding to the 
system's band structure. For systems with a narrow 3d or 4/ band crossing the Fermi level, 
typical of strongly correlated materials, we may divide the polarization into P = P^i + Pr, in 
which Pd includes merely the transitions within the narrow band {3d-3d or ^f-^f transitions) 
and Pr be the rest of the polarization, which includes transitions from the 3d band to the rest 
of the bands and vice-versa. It was noticed that the following quantity can be interpreted 
as the effective interaction among electrons living in the narrow band (Hubbard U): 

U{UJ) = [1 - vPr{i0)]-'^V (3) 

where U can be related to the fully screened interaction W by the following identity: 

W = [l-UPd\-^U. (4) 

This identity explicitly shows that the interaction between the 3d or 4/ electrons is given by a 
frequency- dependent interaction U. Thus the remaining screening channels in the Hubbard 
model associated with the localized d electrons, represented by the d-d polarization Pd, 
further screen U to give the fully screened interaction W. We refer the method of calculating 
the Hubbard U according to ([3]) as cRPA because we have constrained the polarization to 
exclude transitions within the narrow band {d-d transitions). Although the formula in ([3]) 
has been obtained within the RPA, the result is actually exact provided Pr is exact, as was 
shown recently [8|. 

In the following, we retain only the local components of the effective interaction on the 
same atomic site by taking the following matrix element: 

Ul^lm, = J d'rd'r' 02^(r)0z..(r)t/(r, r')0l3(r')0L,(r') (5) 

where cht is a C LMTO |9| orbital {3d or 4/) centered on an atomic site and the interaction 
U{r,r') is the static {u = 0) value of Eq. ([3]). In calculating U we have approximated 0^ 
by the "head" of the LMTO, i.e., the solution to the Schrodinger equation inside the atomic 



sphere. This is expected to be a reasonable approximation because the ( states are rather 
locahzed. LMTO is just one possible choice for the one-particle orbitals but other choices 
are perfectly legitimate. For example, the newly developed NMTO (where N is the number 
of energies chosen to span the region of interest) [10] and the recently proposed maximally 
localized Wannier orbitals |JJj] are possible choices. It is worth noting that the U entering 
the Hubbard model will inevitably depend on the choice of the one-particle basis 0^ defining 
the annihilation and creation operators, no matter what method we use to calculate U{r, r'), 
which is independent of the basis functions used in the band structure method. 



B. LDA-KU 

In the spirit of the LDA+f/ approach we introduce an orbital-dependent exchange- 
correlation operator 

K = XI \(I^RLa)V;^L,R'L'{(pR'L'a\ 
RL,R'L' 

acting among a localized set of electrons. The LMTO head is in general denoted by 
site index R, angular quantum number L = {Im) and spin a. In addition to the usual 



single-particle LDA 



TB-representation 121] we get 



iamiltonian, we include appropriate matrix-elements of V^. In the 
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with 

(0RLcr|X^'L'a) = ^RR'^LV + O Rih^^L ,R' L' ■ 

We have used V^^ /j/^/ = V§^j^5rri5ll' , an assumption which is confirmed numerically, 
ther, the diagonal overlap matrix o as well as the hamiltonian matrix h are given in 12]. 
Consider next V^i^. Assuming a spin- independent Hubbard U , and a diagonal spin-density 
matrix n^^^L' — '^'rl^ll'[^ we obtain: 

L'a' U 

= ULL,L'L'n]il, + {UlL,UL' - ULL',L'L)nRL>. (6) 
L' 

Now Ull^l'l' is substantial for all LL' in contrast to Ull'^l'l which is rather small, except 
when L = L'. We shall use Ull,l'l' = U independent of L,L' and Ull'^l'l = J for L L' 



which result in the simple form 

where the double counting term suggested in Ref. [l| has been added. For a fixed value of 
U and J, the matrix-elements are evaluated and added to the LMTO-Hamiltonian prior to 
diagonalization. The density-matrix ra^^ is updated every iteration using the eigenvectors 
as well as the overlap matrix. The corresponding term, which has to be added to the total 
energy functional, is given by 

RLa 

where = YIrlct'^rl- should be noted that already the simple form of the non-local 
potential gives rise to upper and lower Hubbard bands with an energy separation given by 

iu-J). 

C. Selfconsistent LDA-^U 

The cRPA method requires as input eigenf unctions and eigenvalues (fixed during the cal- 
culation) and delivers as output the Hubbard U matrix. On the other hand, the LDA+t/ 
method needs a [/-matrix (fixed during the calculation) as input and gives as output eigen- 
functions and eigenvalues. The main point of the present work is to merge these two schemes 
in a selfconsistent way. 

We summarize the iterative steps: 

1. Firstly we do a normal cRPA calculation j3] in order to achieve the initial Hubbard 
U matrix (iteration one; matrix Ui) to be used in the LDA+f/ calculation. 

2. After the LDA+f/ calculation has converged we save the output LDA+f/ eigen- 
functions and eigenvalues and use these to calculate U within cRPA (Eqs. ( 12|3|5|) ). in order 
to find the updated f/-matrix for the next LDA+f/ calculation (iteration two; matrix f/2). 

3. The procedure is continued until the f/-matrix is stable i.e after n iterations we 
have Un+i ^ f/„. 
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The size of the [/-matrix is rather large, however many elements are related by symmetry. 



III. RESULTS AND DISCUSSION 

All the results presented in this paper used the simple form of the non-local potential, 
because a substantial number of tests have shown that more elaborate forms of the potential 
do not influence the final results. The most important finding in the present work is indeed 
the possibility to converge the [/-matrix within the defined self-consistency cycle. In all 
cases studied convergency is reached within a reasonable number of iterations. 

To illustrate the applicability of the present scheme to real materials we have applied the 
scheme to NiO, which is an epitome of the charge transfer insulators, and Gd. These two 
systems have been extensively studied both experimentally and theoretically. The NiO LDA 
band gap is known to be too small and likewise the LDA exchange splitting in 4/ Gd is too 
small. These provide a motivation for improving upon the LDA. 

A summary of some results for our prototype systems: For NiO, the self-consistent de- 
termined values U = 6.6 eV and J = 0.9 eV, improves the bandgap (2.5 eV) , compared 
with conventional LDA, though too small in comparison with experiment (4 eV). The ex- 
change spin-splitting of the /-bands in Gd are found to compare rather well with experiment 
(~ 12 — 13 eV) using our self-consistent determined values of U = 12 A eV and J = 1.0 eV. 
We have also calculated the Gd (NiO) magnetic moment to be /i = 7.8 (1.5), which is compa- 
rable to the experimental value /i = 7.6 /ififl^ (1.6-1.9 /i^), and an improvement compared 
to LDA. 

We first discuss Gd, where the LDA+f/ bandstructure corresponding to the self-consistent 
values of U and J are displayed in Figs. ([1112]). The majority (spin up) / bands are centred 
around -11 eV and the minority ones around 3 eV. The occupied spin up bands are very 
narrow due to shielding by the 5s and 5p electrons, due to the hybridization with other 
bands the unoccupied minority bands display some dispersion, making it difficult to extract 
the exchange splitting. However, we estimate that our calculated exchange splitting at 
convergency is somewhat too large by say ~ 1 — 2 eV. 

We note that our parameters differs significantly from those previously used in literature. 



Harmon et al 



15| found U = 6.7 eV and J = 0.7 eV using a supercell approach. The 



experimental gap (splitting between the PES and BIS main-peaks) is given by Eg = Ejy^i + 
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En-1 — '^Ecsi which from purely atomic considerations is predicted to be f/ + 6J, using 
N = 7 spin-up electrons in the groundstate (GS). With the parameters of Harmon et a/. |l5|. 
an underestimation is obtained, resulting in a splitting of 11 eV. As can be seen in Fig. 2, 
the 4f states no longer form a narrow atomic-like band but hybridized with other states in 
the same energy range. Thus, the atomic picture used to estimate the exchange splitting 
may not be valid anymore. 



The frequency-dependent U 



19j from the normal cRPA calculation, i.e., from calculation 



starting from the LDA bandstructure, is shown in Fig. (j3]). We note the dramatic change 
in U for small energies, shooting up to the self-consistent value of U already within 2 eV. 
In fact, the frequency dependence would have become even stronger if we had not used a 
life-time broadening when calculating the response function. Using a tetrahedral method 
for the Brillouin zone integration without a life-time broadening would probably result in 
a decrease in U from its zero-energy value before it shoots up to a large value at around 

1.5 eV. This behavior is in contrast to the transition-metals studied earlier Q]. Towards 
self-consistency we noticed a significant change in U already in the second cRPA calculation 
(1 iteration) ; [/ is in fact enhanced for small energies giving rise to a quite smooth curve 
with weak dependency on frequency. As seen in Fig. ([3]), the frequency dependence of U 
is indeed much weaker after self-consistency, with a relatively constant value of [/ = 12 
eV in the frequency range around 5 eV. The weakening of the energy dependence of U for 
small energies may be explained by the increase in the exchange splitting of the up and 
down 4/ states. As the occupied 4/ states are pushed down the excitation energies from the 
occupied 4/ states to unoccupied states increase. Similarly, as the unoccupied 4/ states are 
pushed up, the excitation energies from occupied states to the unoccupied 4/ states increase. 
Thus, the peak structure in the imaginary part of the screened interaction arising from these 
excitations is shifted to higher energy. Through the Kramers-Kronig relation this results in 
much smoother behavior of U at low energy. This result is very encouraging since it gives 
justification for using a static value of U. 

Finally we consider NiO, where the LDA+U bandstructure corresponding to the self- 
consistent values of U and J are shown in Fig. ( H]). For this system cLDA calculations 
yields U = 8 eV and J = 1 eV pj], which is comparable to our self-consistent values of U = 

6.6 eV and J = 0.9 eV. The gap obtained using the cLDA parameters is 3 eV l|, compared 
to the experimenal gap of 4 eV[16|. The difference between our and the cLDA U (1.4 eV) is 



8 



reflected in our decreased bandgap of 2.5 eV. 



TABLE I: A summary of results for U and J obtained with tlie present method in comparison with 
other methods (in brackets). We compare also the magnetic moments with experimental findings 
(in brackets). 

U (eV) J (eV) Magnetic moment (/xb) 

NiO 6.6 (8.0 [1]) 0.9 (1.0 fl]) 1.5 (1.6-1.9 [17, 18]) 
Gd 12.4 (6.7 [15]) 1.0 (0.7 [15]) 7.8 (7.6 [14]) 



As in the case of Gd, the Hubbard U a.s a. function of frequency undergoes a significant 
change as self-consistency is achieved. Starting from the LDA bandstructure, the resulting 
U calculated using the cRPA method exhibits a strong energy dependence at low energy. 
As the band gap increases, the energy dependence of U at low energy becomes smoother. 
The explanation of this behavior is similar to the case of Gd, namely, as the gap increases 
the peak structure in the imaginary part of the screened interaction is shifted to higher 
energy, and through the Kramers-Kronig relation, it results in a smooth behavior of U at 
low energy. 

The too large 4/ separation obtained in the self-consistent LDA+f/ scheme may arise 
from a shortcoming of the LDA+ U scheme itself rather than in the RPA used in calculating 
U. A similar problem is also observed in the so-called Quasiparticle Self-consistent GW 



(QSGW) scheme [20|. As the 4/ separation becomes larger, the screening associated with 
the 4/ bands becomes weaker, which leads to a larger U. This in turns tends to decrease 
the screening strength and so forth. This indicates a shortcoming of the theory, namely, the 
absence of energy- dependent self-energy and the vertex correction, the latter is also left out 
in the QSGW schemes. 

A further problem that plagues the LDA+t/ scheme is the double-counting problem. This 
problem becomes apparent when the relative position of the correlated bands with respect 



to other bands is important. This relative position is rather sensitive to the double- counting 
formula used in the scheme. We believe this double-counting problem is responsible for the 
incorrect positioning of the 4/ bands in Ce as well as the 3d bands hi^NiO, giving a too small 
band gap in the latter, which has also been found in other works [21]. While the separation 
between the unoccupied Cg and occupied t2g bands of nickel is reasonably well reproduced, 
the relative position of these 3d bands with respect to the oxygen 2p bands is presumably 
incorrect. In the case of the transition metal oxides, such as NiO, the band gap is formed 
between the unoccopied Ni eg band and the occupied O 2p band. 

IV. CONCLUSION 

We have developed a new self-consistent LDA+f/ scheme, in which the important pa- 
rameter U is determined self-consistency using the cRPA method. As test cases we have 
considered NiO and Gd and it is shown that the scheme does yield converged results. The 
exchange splitting in Gd has been found to be too large by 1-2 eV whereas the band gap 
in NiO has been found to be too small, 2.5 eV compared with the experimental value of 
about 4.0 eV. An interesting finding is that the energy dependence of U at low energy is 
found to be much smoother after self-consistency compared with the result obtained from 
the LDA bandstructure. This provides justification for using a static value of U. Our results 
indicate some short-comings of the LDA+ U scheme, in particular the incorrect positioning 
of the 4f states in Gd and the 3d states in NiO points to a problem with the double- counting 
term. Investigating different forms of the double-counting term within the newly developed 
self-consistent LDA+ U scheme could be a fruitful direction to pursue in the future. 
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FIG. 1: Gadolinium spin up bands using the self- consistent determined parameters: U = 12.4 eV 
and J = 1.0 eV. 

Fermi energy at eV and the directions displayed are 1/2(1,1,1) — > P ^ (1, 0, 0). The corresponding 
total DOS and / partial DOS are also displayed. 
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FIG. 2: Gadolinium spin down bands using the self- consistent determined parameters: U = 12.4 
eV and J = 1.0 eV. 

Fermi energy at eV and the directions displayed are 1/2(1,1,1) — > P ^ (1, 0, 0). The corresponding 
total DOS and / partial DOS are also displayed. 



14 




FIG. 3: Frequency-dependent U of gadolinium. 
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FIG. 4: NiO bands using the self- consistent determined parameters: U = 6.6 eV and J = 0.9 eV. 
Fermi energy at eV and the directions displayed are 1/2(1,1,-1) — )• F — )• 1/4(1, 1, 1). 
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